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Abstract 

We propose the autocorrelator of conductance peak heights as a signature of 
the underlying chaotic dynamics in quantum dots in the Coulomb blockade 
regime. This correlation function is directly accessible to experiments and 
its decay width contains interesting information about the underlying elec- 
tron dynamics. Analytical results are derived in the framework of random 
matrix theory in the regime of broken time-reversal symmetry. The final ex- 
pression, upon rescaling, becomes independent of the details of the system. 
For the situation when the external parameter is a variable magnetic field, 
the system-dependent, nonuniversal field scaling factor is obtained by a semi- 
classical approach. The validity of our findings is confirmed by a comparison 
with results of an exact numerical diagonalization of the conformal billiard 
threaded by a magnetic flux line. 

PACS numbers: 73.20.Dx, 05.45. +b, 72.20.My 
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I. INTRODUCTION 



During the last five years the fabrication of ballistic two-dimensional electron gas mi- 
crostructures in GaAs/Gai^ALjAs heterostructures has opened up a new field: the study of 
quantum manifestations of classical chaos in condensed matter systems.0 The first theoreti- 
cal discussion! of the so-called quantum chaos in microstructures dealt with open systems, 
i.e., cavities on the micrometer scale or smaller connected to the external world by leads that 
admit a continuous flow of electrons. Several successful experiments followed. i~Q Nowadays 
the theoretical study of quantum chaos is also well developed for closed microstructures, 
known as quantum dots.! These are cavities coupled to the external leads by tunnel barriers 
that prevent a continuous flow of electrons and make the electric charge inside the dot quan- 
tized. The striking characteristic of quantum dots is the appearance of very sharp peaks 
in the conductance as a function of the gate voltage. These peaks are roughly equidistant 
due to the strong charging effect, but their heights vary randomly by an order of magnitude 
or more. The general interest in this phenomenon, combined with the existence of interest- 
ing experimental data on quantum chaos for closed systems! motivated several theoretical 
studies. These works proposed universal conductance peak distributions for quantum 
dots in the Coulomb blockade regime as fingerprints of the underlying chaotic dynamics. 
Similar predictions for open systems00 have been difficult to observe experimentally^ be- 
cause of the strong influence of dephasing.0 

The initial experiments on Coulomb blockade peak heights in quantum dots! have 
paved the way for the first convincing measurements of statistical distributions of peak 
heights. However, the full distribution of the conductance peaks is not a trivial quan- 
tity to measure. Since each peak sequence is not very long, good statistics demands great 
experimental effort. One possible solution is a true ensemble averaging, which requires many 
different samples. In Ref. [12] it was suggested that one could obtain such an average from a 
single device by changing its shape in situ using external electrostatic potentials. This diffi- 
cult procedure has now been tested experimentally in open@ and closed systems.0 Another 
seemingly simpler solution is to generate more statistics by varying an applied magnetic field 
(taking care to measure at 5-field steps larger than the conductance autocorrelation width 

In order to test more thoroughly quantum manifestations of classical chaos in quantum 
dots in the tunneling regime we propose a different quantity to be investigated: the para- 
metric conductance peak autocorrelation function. Because it is not difficult to cover a 
parametric range which extends over several autocorrelation lengths, correlation functions 
need a small number of peaks and should be easy to obtain experimentally. Moreover, 
correlators of a fixed eigenstate are also of theoretical interest. Analytical calculations us- 
ing random matrix theory beyond a perturbative approach seem to be very challenging. 
In addition, the nonuniversal scaling contains interesting information about the underlying 
classical dynamics. To put the parametric dependence in contact with physical quantities 
and develop a satisfactory theoretical treatment of nonuniversal scaling factors, one has to 
invoke semiclassical arguments. In this sense the study of the autocorrelator of conductance 
peaks is an example where semiclassical and random matrix theories complement each other 
in the understanding of mesoscopic phenomena. 

The first step towards the description of the main features of chaos in quantum dots 
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is to consider the problem of electrons moving in a cavity. We assume that the Coulomb 
interaction within the cavity can be taken into account in a self-consistent way, yielding 
noninteracting quasiparticles. As a result of extensive studies of two-dimensional (2D) dy- 
namical systems, one has learned that when the classical motion in the cavity is hyperbolic, 
the quantum spectrum and the wave functions exhibit certain universal features.^ Through- 
out this paper we shall consider a statistical approach which is valid for chaotic systems and 
also works fairly well for systems whose phase space is predominantly chaotic.@ 

The paper is organized as follows. In Sec. [lj we present our analytical and numerical 
treatment based on random matrix theory (RMT). We concentrate on the case of broken 
time-reversal symmetry (TR), but also show numerical results for the TR preserved case 
(spin-orbit interactions are always neglected). In Sec. |T| we study a dynamical model to 
illustrate our findings. There we conjecture that the semiclassical periodic orbit theory can 
be used to estimate the typical magnetic field scale of any correlator of a quantum dot in a 
realistic situation. We conclude in Sec. [IV] with a discussion on how the presented results are 
robust with respect to dephasing effects and on possible experimental realizations. Several 
more technical details are left to the Appendices. 



II. STOCHASTIC APPROACH 

One of the basic tools to investigate the conductance in mesoscopic devices is the 
Landauer-Biittiker formula.0'@ For a two-lead geometry, where leads are denoted by L 
(left) and R (right), the conductance G is given by 

G = ^E|^(^,0)f , (2.1) 

ab 

where the sum runs over the channels {a} on the lead L and the channels {b} on the lead 
R. The scattering matrix S(E,<f)) connects right and left channels and is evaluated at the 
energy E and magnetic flux <p. The factor of 2 is due to spin degeneracy. 

In the case of quantum dots in the tunneling regime this formula becomes particularly 
simple, since the S matrix can be written as 

S^E, 0) = S° ab (E, 0) - i £ 7w 7 to + 0«r)/A), (2.2) 

v hi — e v + il „/ z 

where the matrix elements 7 0I/ give the probability amplitude of a state in channel a to 
couple to the resonance v in the dot. The total decay width is given by T u = J2 C \lcu\ 2 , 
where the sum is taken over all open channels c. We shall assume that the contribution 
of direct processes is very small and neglect the off- resonance term S°. Equation Q2.2Q is a 
very good approximation to the S matrix when the average decay width (T) is much smaller 
than the mean spacing between resonant states A. This is indeed the case for the isolated 
resonances observed in quantum dot experiments. (Hereafter we use A to denote the level 
spacing caused solely by single-particle states within the cavity, as opposed to the spacing 
between conductance peaks observed in the experiments.) When the condition (r)/A <C 1 
is not met, unitarity corrections to Eq. ( |2.2|) become important and they demand a different 
parametrization of the S matrix (see, for example, Ref. pH) . 
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The matrix elements in Eq. ( |2.2|) can be obtained by invoking the i?-matrix formalism.!^ 
Such an approach was originally proposed in the study of nuclear resonance scattering and 
has more recently been successfully applied in the context of mesoscopic physics.00 The 
partial decay amplitude 7^ is given by 



^Jds X * c (r)Mr), (2.3) 



7c 



where ip u is the eigenfunction of the resonance v inside the dot with the appropriate boundary 
condition and Xc is the wave function in the channel c at energy e v . The integration is done 
over the contact region between lead and cavity A useful quantity for the forthcoming 
discussion is the partial decay width, defined as Y cv = |7a/| 2 - We remark that Eq. (|2.3|), 
as it stands, does not contain barrier penetration factors. Throughout this paper we will 
assume that penetration factors are channel independent and therefore will influence only the 
average decay width (r), but not its distribution. This approximation is certainly valid in 
the case of a single open channel in each lead, the case for which we specialize our analytical 
and numerical results. When the constrictions connecting leads to the cavity are smaller 
than one electron wavelength but the cross section in the leads is not, strong correlations 
among partial decay amplitudes at a fixed resonance exist .Sill For this general situation 
it is still possible to incorporate different penetration factors and the effect of correlations 
into our results, although we do not believe that the final expressions could be cast in a 
simple analytical form. Nevertheless, one would need to obtain all parameters either from 
the experiment, or from a satisfactory dynamical model that comprises both cavity and 
leads. 

Finite-temperature averaging caused by the rounding of the Fermi distribution can be 
done straightforwardly. In the limit T cv <^_kT < A the conductance peak corresponding to 
an on-resonance measurement is given bj 
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G v =^-(^-]g v with g u = VluTru , (2.4) 

where Tl(r)u is the decay width for the resonance v into open channels in the L(R) lead. 
In other words, Tl(r) v = J2 c &l(r)^cu- Equation ( |2.4| ) does not take into account thermal 
activation effects. However, one should have in mind that in the tunneling regime very 
small values of the partial decay widths may occur relatively close in energy to large ones. 
In addition to that, the single-particle spacing fluctuates and can reach values smaller than 
A. As a result, even at low temperatures (less than 100 mK) a very short peak may be 
substantially influenced by neighboring resonances. This effect causes the conductance peak 
to cross over to a different regime, where the height becomes proportional to T rather than 
1/T (even though the inequality <C kT <C A is still approximately satisfied) .0 We will 
avoid further consideration of this or any other similar effect by restricting our analysis to 
peaks whose heights have exclusively the characteristic 1/T behavior. 

In what follows we are going to evaluate the autocorrelator of conductance peak heights 
in terms of an external parameter X which acts on the system Hamiltonian. (At this point 
we do not need to specify X, which could be, for instance, a shape variable or a magnetic 
flux.) In order to obtain the autocorrelator we use a statistical approach that shall be 
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described in detail in the next subsection. Unfortunately, a complete analytical solution for 
this problem is not yet available and presently seems to be a formidable task. Therefore our 
analytical treatment is based on a perturbative expansion for small values of X. We resort 
to numerical diagonalizations of large random matrices to conjecture the complete form of 
the autocorrelator for all other values of X. 

We also derive an exact analytical expression for the joint probability distribution of 
conductance peak heights and their first parametric derivative. This distribution, besides 
being of theoretical interest by itself, can also be employed alternatively in the calculation 
of the small- X asymptotic limit of the conductance peak correlator. 

Before proceeding, we first present some basic quantities entering the calculation. Since 
we are only considering two-point functions, the knowledge of the partial decay widths at 
two parameter space points, X\ and X 2 , is required. We will need to specify r cj ,(X 2 ) only 
up to first-order terms in an expansion around T cu (Xi), namely, 

T CU (X 2 ) = V CV {X 1 ) + Xk^Xt) + 0(X 2 ), (2.5) 

where X = X 2 — Xi . The coefficient A is given by first-order perturbation theory in terms 
of eigenfunctions, eigenvalues, and matrix elements of the perturbation. Assuming that the 
Hamiltonian may be written in the form H(X) = Hq + XU and using Eq. ( |2.3| ), we find 
that 

Ao/ = 2^ —— » ( 2 - 6 ) 

with i? ^v(r) = e v ip u (r) and all partial decay amplitudes, eigenvalues, and matrix elements 
evaluated at X\. So far no statistical assumption has been made. 



A. Parametric conductance peaks autocorrelation function 

The starting point of our discussion of the conductance peak autocorrelation function is 
Eq. (|2.4j ). In this subsection we shall mainly treat the situation of TR broken by an external 
magnetic field. We identify the magnetic flux with the previously introduced variable X. 
In the spirit of RMT,lll the analysis that follows will lead to a universal curve for the 
conductance peak autocorrelation, which is expressed in terms of a scaled X and the coupling 
of the dot to the external leads, denoted by the parameters a^L- Both a^L and the typical 
scale of X are system-specific quantities, depending on the dot shape, the quality of the dot- 
lead couplings, and the Fermi energy. Ideally, the scaling of X should be extracted from the 
experimental resonance (level) velocity correlator.il When this information is not available, 
it is not straightforward to put X in correspondence with experimental parameters. To 
overcome this problem we invoke a dynamical system and shall discuss it at length in Sec. 
III. Similarly, to obtain or^ one needs to model specific features of the tunneling barriers. 
We avoid this procedure by expressing our results in terms of the mean decay width (r), 
which then becomes a fitting parameter in comparisons with experiments. 

The autocorrelator of peak heights is defined as 

C„( X) s <* (X - f ) *(* + f )) - (* (X - f ) ) <* (X + f )) , (2.7) 
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where (• • •) denotes averages over the resonances v and over different values of X, but will 
later be interpreted as an average over the Gaussian unitary ensemble (GUE) of random 
Hamiltonians.EJ The difficulties involved in the exact calculation of C g (X) become clear once 
we expand the conductance peak height in powers of X and assume translation invariance 
(independence of X), 



N ( 1^\ n X^ n 



c 9 (x) = E 



n=0 



d n g„{X) 
dX n 



{g v {X))\ (2.8) 



(Notice that odd powers disappear since the autocorrelator is by construction an even func- 
tion of X.) As presented in Eq. (|2.8|) , the task of finding the complete functional form of 
C g (X) requires not just the knowledge of the second moment of all derivatives of g v (X), but 
also the summation of the Taylor series. Presently we do not know of any other alternative 
method that could lead to the complete analytical form of C g (X). An analogous situation is 
encountered in the calculation of the more commonly studied level velocity autocorrelator^ 2 ] 



(see the detailed discussion of Ref. |33|). The small X asymptotic limit can, however, be 
evaluated exactly through standard random matrix methods. Moreover, an approximate 
curve for all values of X can be found through numerical simulations. 

We start with the small- AT asymptotics. For this purpose we keep in Eq. fl2.8|) only the 
two lowest order terms, 

C g {X) = (7,(0) + C"J0)¥- + 0(A 4 ) , (2.9) 



with 
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C,(0) = (( ,:"Tr' ) ( , ' /,Vl ) (2-10) 




and 



C"(q) = l( T Rv T Lu X 2 /A Ru ^ A Lu A Ru + A Lu \ 2 \ (2 11) 

\ V + r^j, / vr^, Yl v t Rv + Tl u ) j 

Now we introduce a statistical model: we assume that due to shape irregularities and the 
presence of a magnetic field, the Hamiltonian can be modeled as a member of the GUE.0 
Under this assumption, C g (0) has already been obtained. 00^1 For example, in the sim- 
plest case, where one has two equivalent uncorrelated open channels (one in each lead), 
C g (0) = 4(r) 2 /45. The second-order coefficient requires a more complex treatment, since 
Eq. Q2.ll ) involves not simply eigenfunctions but eigenvalues as well. For Gaussian ensem- 



bles, a great simplification is possible due to the statistical independence of matrix elements 
of the Hamiltonian. This will allow us to divide the ensemble average into averages over 
eigenfunctions and eigenvalues. 

For the sake of simplicity, we will take the matrix elements of U as Gaussian distributed,^ 
with 



a 2 



(U^)=0 and (U^U Xp ) = , (2.12) 
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where a is a measure of the perturbation strength and N is the number of eigenstates 
(resonances) considered. It immediately follows that 



2^ 

N 



\lcfj,\ \lcv\ 



(2.13) 



c£L(R) n^u 



'Ufl 



and 



(A 



Ru^Lu)u = ~jy E 

cdR,c'aL fij^u 



E 



1 'cvT / C 'vl 'cfil ' c ' '/, 



+ C.C. 



(2.14) 



'UfJ, 



with e Vfl = e v — e^. Up to this point the formulation of the problem is very generic and no 
restriction has been made either on the number of open channels or on their transmission 
coefficients. A quick inspection of Eqs. ( |2.13j ) and ( |2.14j ) indicates that the solution of 
the most generic case (multichannel leads) is possible, but the algebra is very tedious and 
intricate. Recent experimental realizations of quantum dots in the tunneling regime@Jii 
indicate that one typically has only one open channel in L and R. Owing to that, we limit 
our calculations to the one-channel case. 

Carrying out the average over U by using Eqs. ( |2.13| ) and (|2.14j ) to simplify Eq. Q2.11| ), 
we get 



N 



Ru 1 - Lu 



Fru + ^Lu 



1 

F 2 " 

^U ^UfM 



E 



1 Rfi 

r 3 

1 Ru 



+ 



L/i 



r 3 



+ 



Lu 



iRvlLulRylLp 

r 2 r 2 

1 Ru L Lu 



+ C.C. 



(2.15) 



As pointed out before, the independence of eigenvalues from eigenfunctions allows the fac- 
torization of the right-hand side of Eq. (|2.15| ) into two decoupled averages, namely, 



2^ 
N 



A T A £ 



(2.16) 



The first factor is given by 
A r = ( 



T R ,T 



R,u L Lu 



Ru 



Lu 



f 3 " 
1 Ru 



r 3 

1 Lu 



(2.17) 



{r R ,r L } 



where the average {T} is taken over both leads L and R and resonances v and /x with v ^ fi 
(terms with an odd number of distinct amplitudes drop out). The second factor is 



A, 



(2.18) 



and it involves only an average over the eigenvalues. Since the lead contacts are typically 
farther than one wavelength A apart, the decay widths at R and L are assumed to fluctuate 
independently. Hence, the calculation of A r requires the convolution of four \ 2 distributions 
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with two degrees of freedom,!! P2CO = exp(— Y / (r)) (notice that we are considering a single 
open channel per contact and identical leads), and we obtain A? = 2(r) 2 /15. 

The eigenvalue average in Eq. ( |2.18| ) is more complicated to carry out and needs some 
additional assumption. Because we are going to take the thermodynamic limit (N — > 00) 
later on, effects due to variations in the density of states are negligible. Consequently, we 
may place the reference eigenstate e v at the center of the spectrum and write 



{<*} 



where p(E) = (£^ =1 8{e v - E)) is the average density of states. Before proceeding to 
evaluate this quantity in the RMT framework, let us first discuss an illustrative and extreme 
situation, the picket fence spectrum with spacing A. In this limit the averaged sum simplifies 
to (2/A 2 ) Y^=\ n ~ 2 , which readily yields A e = 7r 2 /3A 2 . Because the Gaussian ensembles 
hypothesis implies fluctuations in the spectrum, the actual value of A e for the GUE should 
be larger but of the same order of magnitude. 

The average over eigenvalues in Eq. (|2.19|) can be restated as an average over a spectral 
determinant, 



If ( N N \ N 

Ae = (c\\7 d£ i--- de N exp -— el + 2 £ In \e V(t \ £ J2 



5(£u) 



det(H')tr(H- 2 )) u , (2.20) 



p(0) Z N \ N J \ v ' v 'Is. 

where H is an (N — 1) x (N — 1) GUE matrix and the normalization constant is given byll 



r I N N \ 

J ds\ • • • d,£pf exp -^X/ £ " + 2 zZ ln l e ^l 



/\2\N 2 /2 N 

= &*)W\±) Uk\. (2.21) 

[Notice that A = l/p(0) = X/ttN is the mean level spacing at the center of the spectrum.] 
In the limit of N ^> 1 Eq. ( |2.2C| ) can be evaluated by the fermionic methocfl (see Appendix 
^ for details) and one finds that 

2tt 2 

A = ^ ■ (2.22) 

e 3A 2 V ) 

As expected, the actual value of A e is larger, but still in fair agreement with the picket fence 
estimate. 

The next step is to rescale the perturbation X to a dimensionless form, in such a way 
that all system-dependent parameters are eliminated. One possibility is to use {(dg/dX) 2 ) 
as the rescaling parameter. As will become clear later on, we do not find this procedure very 
interesting from the physics viewpoint because this quantity cannot be easily calculated 
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given the underlying dynamical system. We rather follow an idea originally proposed by 
Szafer, Simons, and Altshuler: Recall that the perturbation strength a (the nonuniversal 
scale in the above calculations) also appears in the level velocity correlator C^(X)^fi 

when evaluated at X = 0, namely, 

C M = = . (2 .24) 



The statement implicit in the original works of Refs. |32| and |33| is that the quantity y C„(0) 
sets the scale for any averaged parametric function (f(X)), provided that the system dy- 
namics is chaotic in the classical limit. In Sec. |ITT]we will show that C v (0) can be obtained 
by semiclassical arguments, once details of the confining geometry of the dot are known. 
Therefore, in analogy to their analysis of level velocity correlator, we apply the rescalings 



x = XyC v (0) and c g (x) = C g (X)/C g (Q) to arrive at the following universal (dimensionless) 
form: 

Cg(x) = 1 - ttV + 0(x 4 ) , (2.25) 

valid for x C 1. 

To extract the nonperturbative part of well as its x ^> 1 asymptotic limit, 

we relied on a numerical simulation. We performed a series of exact diagonalizations of 
random matrices of the form H(X) = Hi cos(X) + H 2 sm(X), with Hi and H 2 denoting 
two 500 x 500 matrices drawn from the GUE. This model for the parametric dependence is 
rather convenient for the simulations and later data analysis because it does not make the 
level density depend on X, nor does it cause the eigenvalues of H(X) to drift with X. It 
is helpful to think of [H(X)]ki as the matrix element of the Hamiltonian in a discrete space 
representation. As a result, for the one-channel lead case we may simply equal decay widths 
to the renormalized ((r) = 1) wave function intensities at a given point, 

T ku = N\Mk)\ 2 , (2.26) 

where k is the site number and v the eigenstate label. In this way, we were able to generate 
more than 10 5 different configurations of the two-lead geometry, out of which only 2 x 10 4 
were used (we stress that wave functions taken at different v or k are statistically independent 
when the size of the matrix is large enough) . 

For each realization of Hi and H 2 we varied X in the interval [0, 7r/8] and considered 
only the 100 central eigenstates in order to avoid having to unfold the spectrum. In total 
we have run 50 realizations. The final result is presented in Fig. |I|. For comparison, we have 
also shown in the inset the result obtained for the decay width correlator, 

Cr(x) = (r kv (x - 1) r ku (* + y) ) - (r to (x - 1 ) } (v ku (x + 1 ) \ . (2.27) 
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We remark that this correlator is not directly accessible to experiments in quantum dots 
(recall that conductance peak widths are dominated by the thermal rounding of the Fermi 
surface) .il Here we have introduced it with the unique purpose of checking the reliability 
of the numerical simulations in the i>1 range. Contrary to the situation for c g (x), both 
i < 1 and x ^> 1 asymptotic limits of the rescaled Cr(x) = Cr(X)/Cr(0) can be evaluated 
analytically. One finds that (see Appendix |B|) 

/ 1 - 2ttV/3 for x < 1 , . 

cr(x) — > < _ xo ' , . 2.28 

K ' \ 1/(ttx) 2 for x > 1 y J 

The data obtained from the random matrix simulations indicate that the asymptotic tail of 
c g (x) is well described by an x~ 2 law. In the light of Eq. ( |2.28p , this is not very surprising: 



If one lead were more strongly coupled to the cavity than the other, say, r# ^> Tl, we would 
have that c g (x) = cr(x) + 0(Tl/Tr) and therefore c g (x ^> 1) ~ x~ 2 in leading order. We 
cannot rigorously prove, though, that this asymptotic form is also exact when right and left 
leads are identical. 

To conclude this subsection, we briefly discuss the universality class of preserved TR 
without spin-orbit coupling (the case when spin-orbit coupling is present, the symplectic 
ensemble, will not be discussed since it is not relevant to semiconductor heterostructures) . 
The simplest experimental realization is a measurement of the evolution of conductance peak 
heights as a function of shape deformation in the absence of a magnetic field. The general 
approach is the same as above and we assume that the system Hamiltonian can be modeled 
as a member of the Gaussian orthogonal ensemble (GOE).0 However, the i -C 1 asymptotics 
of the conductance peak correlator is now more difficult to calculate. This is not a daunting 
problem because, as seen above, the numerical results reliably recover the correct behavior 
of the correlator for small values of x. Consequently, for the TR-preserved case we relied 
entirely on numerical simulations and did not attempt any analytical calculation. We used 
the same parametric dependence of H(X), but this time drew Hi and H 2 from the GOE. All 
other steps were identical to the GUE simulation. The resulting correlation functions (after 
the proper rescalings) are shown in Fig. |2|. Notice that the large- x asymptotics of both c g (x) 
and Cr{x) are well described by an x~ 2 law, as for the GUE simulations. Here, analogously 
to the GUE, we do not know how to prove analytically that this decay is rigorously true for 
c g (x); on the other hand, we do know that this power law decay is indeed exact for cp(x)S 
We emphasize that the most important characteristic of the TR-preserved, universal c g (x) 
as compared to the TR-broken one is the larger decay width. The widths at half maximum 
height differ by approximately 20%. 

The power spectra of conductance peak height oscillations [the Fourier transform of 
c g (x)] for both GUE and GOE are shown in Fig. |3|. Notice that the behavior is exponential 
only over a small range of k values. One then concludes that a Lorentzian can only be used 
as a rough approximation to the exact curves. Besides, a Lorentzian cannot accommodate 
simultaneously the small and large asymptotic limits of c g (x) presented above in either 
TR-preserved or TR-broken cases. 
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B. The joint distribution of conductance peak heights and their first parametric 

derivative (unitary ensemble) 



In this subsection we evaluate the joint probability distribution 

dg v (x)\\ 



Q(g,h) = (5(g-g v (x))5[h 



dx 



^2.29) 



where x = ^C v (0)X and the Hamiltonian belongs to the GUE. This distribution, although 
not suitable for direct experimental investigations, allows one to easily evaluate higher mo- 
ments of the conductance peak height and its first parametric derivative; in particular, one 
can use Q(g,h) to obtain the two first coefficients in the expansion of Eq. (|2.9|) , 



0,(0) 



and 



CM 



dg / dh g 2 Q{g, h) 

J — oo 



I poo poo 

—— / dg dh h 2 Q(g, h) . 

„ U JO J-oo 



(2.30) 



0,(0) Jo 

We begin by recalling Eq. ( |2.6|) and introducing the coefficients A u into Eq. (|2.29|) , 

^Rv^L 



2.31) 



Q(g,h) = (5[g 



Rv ' J- Lv , 



\8(h-M u 



where 



C„(0) 



Ar Al Ar + Al 



Taking the Fourier transform of Eq. ( [2.32j ) with respect to h, we get 

FrvTlu 



(2.32) 



(2.33) 



Q(g,t) = {s[g 



Tr„ + T 



Lv 



exp(i£M„ 



(2.34) 



Following the same assumptions of the previous subsection, we break up the ensemble 
average into four partial averages. Recalling Eqs. ( [2.13|) and (|2.14j ), we first carry out the 
average over the perturbation U: 



(exp(itM u )) u = Y[ expl 



gH 2 A 2 



iRvlRn + iRvlLp, 



r 2 

1 Rv 



r 2 

1 Lv 



(2.35) 



H^v \ v^i 

Next, we evaluate the average over the decay widths {7^,7/44} forfi^v only. 
I ' gH 2 A 2 _ Yrv1r» + livlL^ \ \ 



exp 



'V\l 



r 2 

1 Lv 



1+g _w±( i + i 



i -i 



r 3 r 3 

1 _Ri/ 1 Lv , 



[2.36) 
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At this point, instead of also taking the average over the remaining decay widths, we consider 
first the average over the eigenvalues {s u }, v = 1, • • • , N. Here we use the following relation, 
valid in the large- N limit and for e u at the center of the spectrum:!! 



n 



1 + 



'Vfi 



B(k) 



(2.37) 



where B(k) is the Fourier transform of 



Bis) 



35 + 14s 2 + 3s 4 
12tt(1 + s 2 ) 4 



^2.38) 



[The evaluation of the average in Eq. ( p.37|) requires a generalization^ of the approach used 
to derive Eq. fl2.22| ). A brief description is given in Appendix y.] After inserting (|2.37|) into 
Eq. ( |2.36|) and inverse Fourier transforming the result we get 



Q(g,h) = (5ig- 



r B „r 



Ru 1 - Lu 



X B 



( 



-3 



^2.39) 



which is more conveniently expressed in the form 

e -4s/<r> 



Q(g,h) 



--H, 



2n(T^^g^f 9/{T \2n^/(T) 



h/(T) 



[2A0) 



The remaining average over the reduced widths an = r^/(T) and — Tl u /(T) appears 
only in the evaluation of the bell-shaped function 



roo roo 

H p {q) = da R da L e ^ aR+a ^ 
Jo Jo 



Hp- 



p*{a R 3 + a L 3 ) 



Q 



p*(a R 3 + a L 3 ) 



(2.41) 



l R ' ) N V 1 v R 

In fact, one can easily carry out one of the above integrations and arrive at 



roo 

HJq) = 2p du e 
Jo 



u(u+l) V 1 



lu + 4 



u 



[2A2) 



It is possible to represent H p (q) in terms of special functions, but we did not find it partic- 
ularly clarifying and therefore do not do it here. The asymptotic limits follow directly from 
the integral representation shown above. For q = 0, we obtain 



H p (0) 



(35/67r) for p — > 



' \ (70/3)yJp/7r for p 
whereas for a fixed p ^ 0, we have that 



oo 



(2.43) 
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H P (q) ~ 0(g~ 4 ) for g > 1 



(2.44) 



Notice that Eqs. ( |2.40| ) and Q2.42j ) together immediately lead to the same x C 1 asymp- 
totics for the peak height correlator shown in Eq. (|2.25|) . From Eq. (|2.44j ) it is obvious that 
(h n ) diverges for n > 2. This can be ultimately related to the level repulsion present in 
the spectrum: Strong anticrossings of levels can cause anomalously large variations of the 
conductance peak height as a function of the external parameter X. When this happens, 
one finds that h ~ where u = \s v +i — <C 1. Since the probability of this event 

goes as P{oS) ~ uj 2 for the unitary ensemble,^! one finds that Q(g, h) ~ l/h 4 for h 1, in 
agreement with our exact calculation. 



III. DYNAMICAL MODEL 

The aim of this section is to compare the results of the previous section with exact 
numerical diagonalizations of a dynamical model. The essential characteristic of a dynamical 
model for this type of study is a fair resemblance to the actual experimental conditions, 
combined with its adequacy to numerical computations. For this reason we chose the (two- 
dimensional) conformal billiard penetrated by an infinitely thin Aharonov-Bohm flux line 



carrying a flux of 0. This model was originally introduced in Ref. [41] and later adopted in the 
study of statistical features of conductance in quantum dots.Hl Using complex coordinates, 
the shape of the billiard in the w plane is given by \z\ = 1 in the following area preserving 
conformal mapping: 

. , z + bz 2 + ce i5 z 3 , n , . 

w(z) = — , 3.1 

V ; VI + 2b 2 + 3c 2 V ; 

where b, c, and 5 are real parameters chosen in such a manner that |ttf'(z)l > for \z\ < 1. 
The classical dynamics of a particle bouncing inside this billiard is predominantly stochastic 
and is unaffected by the presence of the flux line. To describe the flux line we chose the 
following gauge for the vector potential A: 



A(w) = ^ 9f{w) 
y ' 27r\dIm(wY dRe(w) 



, (3.2) 



where f(w) = In [|z(it?)|]. This particular gauge, combined with Neumann boundary condi- 
tions, permits a separation of the Schrodinger equation into polar coordinates (r, 6) of the 
complex parameter z = re %e . The eigenstates ip u thus obtained correspond to the resonant 
wave function appearing in Eq. (|2.3| ).Ej In our numerical treatment the wave function \c i n 
the lead is equal to a transverse sine function multiplied by a longitudinal plane wave. 

In practical calculations]!! one fixes the value of the flux </> and uses as a truncated 
basis the lowest (in our case 1000) eigenstates of the circular billiard {b = c = 0). These 
have the form Jy^nv^G 116 ', where J v is the Bessel function of fractional order v and 7 ni , is 
the nth root of its derivative, J' v {j n u) = 0. The dependence on the magnetic flux enters 
through v — \l — <f>/<f) \, with O as the flux quantum h/e. To solve the Schrodinger equation 
one has to calculate several thousand matrix elements of the Jacobian J = \w'(z)\ 2 in this 
^-dependent basis. This operation is very time consuming. Changes in shape are not a 
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major obstacle, since b, c, and 5 act as prefactors to the matrix elements and no further 
calculation is necessary. (For instance, in Ref. [12| only two different values of <fi needed to be 
used.) In the present work, however, we wanted to change the flux to simulate the simplest 
experimental setup and this required the use of 76 different values of <fi. The spectrum is 
seen to be not only 2ir periodic in 0, but also symmetric around 4>/4>o = 1/2. For this reason 
and, furthermore, to avoid the special points and 1/2 we let <f)/4>o vary in the interval 
[0.1, 0.4]. To circumvent the overwhelming problem of too large amounts of computing time 
we employed the following strategy. We calculated the Jacobian matrix J for only seven 
different values of a = 4>/4>o, namely, from 0.10 to 0.40 in steps of 0.04. Then, for any other 
value of a, the matrix element Jij{a) was found by polynomial interpolation. This was 
checked to give a relative error of at most 10~ 7 . To improve the statistics we also calculated 
the eigenstates for five different shapes by keeping b = c = 0.2 and letting 5 = kir/6, with 
k = 1,2,3,4,5. The spectra corresponding to these shapes are statistically uncorrelated.il! 
Due to the truncation of the basis only the lowest 300 of the calculated 1000 eigenstates were 
accurate enough to be used in the analysis. We also discarded the lowest 100 eigenstates 
because of their markedly nonuniversal behavior. We should emphasize that it is not a 
trivial task to increase the number of usable states. For the asymmetric conformal billiard 
no symmetry reduction of the resulting eigenvalue problem is possible. The presence of 
the flux line constraints the method of analysis to the diagonalization of large Hermitian 
matrices, limiting the number of eigenstates that can be treated efficiently. 



A. Correlation functions for the billiard 

In the following, we present the numerical results for the correlators of level velocity, 
decay widths, and conductance peak heights. (Recall that at present only the last can be 
directly measured in real experiments.) 

Figure § shows the level velocity correlator defined by Eq. (|2.23|) rescaled according to 
c v (x) = C v (x)/C v (0), with x = a/yCJfi). The plotted data were obtained by averaging 
over 76 equidistant values of a in the interval [0.1, 0.4], over the eigenstates between 201 and 
300, and over the fivedifferent shapes mentioned above. We observe a good agreement with 
the analytical results^ for small values of x and with random matrix simulations in general. 
A thorough discussion of the scaling factor <JC v (0) dependence on energy and billiard shape 
is postponed to the next subsection. 

Next we employed the billiard model to obtain the decay width autocorrelation function 
cr(x), Eq. ( |2.27| ). The decay widths Y ku for the billiard were calculated as described in 
Sec. [TJ. The position of the leads and their widths were specified in the following way. 
Based on autocorrelation calculations of the eigenfunction ip v along the perimeter of the 
billiard, we found that the spatial decorrelation for the relevant levels (levels 201 to 300) 
takes place over a distance around 1/60 of the perimeter. We therefore decided to take 
the width of the leads to be 1/24 of the perimeter, i.e., 2.5 times the decorrelation length, 
yielding 24 adjacent leads. To improve the statistics we used all 24 lead position for each 
ip v . Due to the relatively large width of the leads, adjacent leads are not correlated, as was 
verified by obtaining the same result (with larger fluctuations) using only every second or 
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every third lead. The result of the calculation is shown in Fig. [5|. A fair agreement with 
Eq. ( g2p is noted - 

Finally, we calculated the conductance peak correlation c g (x), Eq. (|2.7| ), for the conformal 
billiard using the above mentioned decay widths. We chose all possible pair configurations 
among the 24 lead positions. The result is shown in Fig. |6] and, again, a fair agreement 
with the predictions of Sec. |II A| is observed. Notice that the data for the conformal billiard 
are not fully consistent with a Lorentzian if we fix the a; < 1 asymptotics of the curve 
to be identical to Eq. ( |2.25| ). A squared Lorentzian does not seem to provide a better 
approximation either, although in a recent experiment^! c g (x) was measured and such a 
curve was fitted to the data. It would be interesting to check how well our result for c g {x) 
based on random matrix calculations matches the available experimental data without any 
fitting parameter (in Sec. |Ti [13] we will present a way to predict the typical field correlation 
scale). 

For all three correlators we have noticed large statistical fluctuations between data taken 
at different billiard shapes. We found that most levels around the 300th one (the upper limit 
of reliability in our calculation) still do not show more than one full oscillation within the 
range of flux allowed by symmetry. The averaging over shape deformation was thus crucial 
to get rid of the remaining nonuniversal features. After averaging over the five values of 5 
mentioned previously (see also the following subsection), we found that C v (0) near the iVth 
level obeys the law C v (0) « 1.202^. 

We ascribe the small mismatch between theory and numerics, particularly at x > 1, to 
poor statistics. As mentioned before, the only way to circumvent this problem is to compute 
higher eigenstates. Fortunately, this difficulty does not appear in real experiments where 
the magnetic flux is extended over the whole area of the cavity and the dependence of is 
not periodic. 



B. Energy dependence of C v (0) for billiards 

One of the important features of a dynamical model is that its quantum fluctuations 
display a marked dependence on energy. A quantitative understanding of the field scale of 
the fluctuations and its dependence on energy is important to put any random matrix result 
in contact with measurements in quantum dots. For this purpose, the semiclassical approach 
can be used in a relatively simple form. 

The aim of this subsection is to discuss the energy dependence of the level velocity 
correlator C v (0), which is a measure of the quantum fluctuations, and show that one can 
successfully estimate this quant ityusing exclusively classical quantities. It was already prop- 
erly noticed by Berry and Robnikcil that the typical flux necessary to induce a crossover from 
GOE- to GUE-like spectral fluctuations in a chaotic billiard threaded by an Aharonov-Bohm 
flux line scales with the energy as E 1 ^. The origin of this dependence has a simple semi- 
classical explanation which was nicely worked out by Ozorio de Almeida and coworkers.l^H 
The nonuniversal scaling C v (0) can be obtained in an analogous manner. 

Our point of departure is a recent work by Berry and KeatingE! based on the Gutzwiller 
trace formula. Most of our semiclassical considerations follow their findings. However, 
our interpretation and the method we use to quantify C v (0) are different. To make the 
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exposition self-contained, we shall briefly describe points of their work which are relevant to 
our discussion and comment when necessary. 

The initial step is to approximate C v (<p) by a two-point correlator. (The nature of the 
approximation is evident, since to track down the parametric evolution of a single eigenvalue 
for any is a task that cannot be exactly achieved by considering only Green's functions 
with a finite number of points.) As in Ref. ^2], we write 

(3.3) 

where N V (E, 0) = J2u ®r)(E — E u ((j))) is a smoothed cumulative level density which counts the 
number of single-particle states up to an energy £ at a magnetic flux 0. For the smoothing 
it is convenient to adopt the form dQ v (E)/dE = lm(E — it])^ 1 /?!, where the parameter r] 
is chosen to be much smaller than A. (We would like to point out a change in notation: 
Throughout this subsection angular brackets (• • •) will denote energy averages in distinction 
to the previous ensemble averages.) The energy average in Eq. ( |3.3j ) is taken over a range 
SE around E. Ordinarily, the average over the magnetic flux is taken over a window <50 
in flux which corresponds to little change in the classical dynamics but is semiclassically 
large, i.e., it corresponds to sizable fluctuations in the spectrum. For billiards threaded by 
Aharonov-Bohm flux lines this is not an actual constraint, since flux variations have no effect 
on classical trajectories. Nonetheless, when considering the correlator of Eq. ( |3.3|) one has, 
in principle, to avoid pertaining to the TR-breaking crossover regime. 

When levels are much farther apart than rj, it is straightforward to show that dN^^) / d<p 
can be approximated by A _1 cfe u (<f))/ 'dip. For larger than a certain C such that the cor- 
relation between level velocities of different states /i ^ v is much weaker than for fi = v, 
F v ((j),E) is equivalent to C v (4>) and independent of 77, provided that the energy levels are 
taken to be within SE. This equivalence also holds in the limit of = when the spectrum 
is nondegenerate. In summary, F v (E,<j)) is a good approximation to the level velocity cor- 
relator C v ((p) only for = and > C . In fact, F v (E,(j)) is also the quantity evaluated 
analytically by Szafer and Altshulerll in the context of disordered metallic rings, using a 
diagrammatic perturbation theory based on impurity averaging. 

Following Ref. |42|, we now turn to the semiclassical treatment. In particular, we specialize 
the results to billiards threaded by a flux line. This simplifies the problem enormously since 
if S n is the action of a periodic orbit n, upon applying a magnetic flux we find that 
S n — > S n + 27r/k(; n 0/0o, where w n is the number of times the orbit n winds around the flux 
line. The cumulative level density is expressed semiclassically using the Gutzwiller trace 
formula and one writes the correlator as 




F,(0, E)=[^) ( £ \A m A n \ <ex P U (^j^ + 2™n^) - v{Tn + Tm] } 5. 



w n w m 

I SE 



(3.4) 

where the amplitudes A n (E) contain information on the stability of the orbit n as well as 
its Maslov index. The smoothing of the staircase function gives rise to the exponential 
damping factor rj times T n (E), the period of the closed orbit n. In the semiclassical limit, 
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N(E) ^> 1, the flux range of the TR-breaking crossover is much smaller than 0o- Therefore, 
without affecting considerably the calculation for the pure TR-broken case, we take the limit 
of 5(f) — > 0o- 

Now one arrives at one of the delicate points of the semiclassical approach. The correlator 
F v (<p,E) is expressed as a sum of diagonal and off-diagonal contributions. In contrast to 
the fact that it is a settled matter how to compute the diagonal part, the evaluation of 
the off-diagonal term is still an unsatisfactorily solved problem. It seems that for F v it is 
reasonablec3 to neglect the off-diagonal contribution and we will do so hereafter. 

The diagonal part of F v (<p, E) in Eq. ( |3.4| ) is still difficult to evaluate since, in principle, 
it requires the knowledge of the full set of periodic orbits up to the cutoff h/r]. This can 
be simplified by the following considerations. For chaotic systems, the number of periodic 
orbits grows exponentially as a function of their length L and ergodicity ensures that as the 
orbits become longer they tend to explore the phase space more uniformly. Thus, one can 
define a critical length L c corresponding to an uniform coverage of the phase space by the 
periodic trajectories. For a fixed energy E, L c determines T c , the time when the Hannay 
and Ozorio de Almeida sum rule^l is applicable, allowing the sum over periodic orbits to be 
calculated as an integral over orbital times. Applying this sum rule to the diagonal part of 
Eq. ( |3.4D one obtains 

^• E) "Jtd ^)co S ( 2 ™ ( r)£) exp ) . (3.5) 

Here the overbar stands for an average over the phase space, or, in practice, over an ensemble 
of trajectories. The replacement of the summation over periodic orbits by an average over 
trajectories makes the problem amenable for a computational evaluation of F^ iag (0, E). Such 
a procedure has already been successfully used to estimate the GOE to GUE transition 
parameter in the stadium billiard.^ 

By writing the winding number as an integral over the angular velocity, w(T) = 
dt6(t)/2ir } it is simple to see that for an ergodic system w{T) = 0. The winding number 
variance is given by 

"cspf*' ^ (3 ' 6) 

where C(t) = 6(t')8(t' + 1). For chaotic systems in general the correlator C(t) decays suf- 
ficiently fast in T to assure the convergence of the integral in Eq. (|3.6| ).ii The knowledge 
of the winding number distribution P(w,T) allowed us to evaluate the phase space average 
in Eq. (|3.5| ). To determine P(w,T) for the conformal billiard we have randomly chosen 10 4 
initial conditions and computed trajectories up to 250 bounces for the particular deforma- 
tion b = c = 0.2 and different <5's. We have generated a histogram recording w and T every 
time a trajectory winds around the flux line located at the originB The results displayed in 
Fig. [7| confirm Eq. ( |3.6|) . The variance of w(T) is better written asEl 

/ op \ 1/2 
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where A is the billiard area and ft is a system-dependent quantity computed for the scaled 
billiard with unit area and trajectories with unit velocity. Moreover, our numerical results 
give us confidence that for any given time P(w,T) is a Gaussian distribution (see Fig. ||). 
This is in agreement with the conjecture of Ref. |4l] for periodic orbits. 

Substituting this result into ( |3.5| ) and baring in mind that in the semiclassical limit 
r]T c /h <C 1, for cf) = we obtain 

F^(0,E)*j 2 . (3.8) 
Using the leading term in the Weyl formula, N(E) ~ AmE/2nh 2 , we can write 

0c = — T77 • (3.9) 

Notice that C <C <p . 

Although F(4>,E) fails as an accurate approximation for C v {4>) over the entire range of 
<p, we expect it to work for = and <p ^ <Pc- Therefore, the semiclassical quantity 4> c given 
above should yield a good approximation to the exact inverse field scale JC v (0). We used 
the same procedure described to obtain P(w, T) to compute k as a function of S. The results 
are show in Fig. |5[ Figure |H] shows C v (0) as a function of iV from quantum mechanical 
(Sec. [Ill A[) and semiclassical (<f)~ 2 ) calculations. From the numerical diagonalizations we 



found that the proportionality factor between C v (0) and N 1 ^ 2 varies between 0.94 (5 = 5tt/6) 
and 1.54 (S = vr/6), while the semiclassical estimate gives 0.90 for 5 = tt/2, for instance. 
Clearly, both quantum and semiclassical calculations indicate that C v (0) depends on the 
billiard shape. When comparing the results of these two calculations, one should also note 



the size of the large fluctuations in the data presented in Fig. |10[ As already remarked in 
the previous subsection, these large statistical fluctuations are due to the limited data set 
used in the simulations. 

To put the semiclassical result in direct contact with experiments, one needs some system- 
specific information to compute k. In the above discussion, k gives a measure of how fast the 
winding number variance increases with time. For the more physical situation of extended 
B fields. /« measures the rate of increase in the variance of accumulated function 
of time.c2l It is interesting to notice that the semiclassical interpretation of the scales of the 
conductance autocorrelation function for quantum dots and open cavities is very similar. 
Since quantum dots have T A, the escape time %/Y is always semiclassically large. 
In this regime, the physics is dominated by the classical decorrelation time r implied in 
Eq. ( |3.6| ). The situation is very different for open cavities, where the escape time plays an 
important role, since it is comparable with r. For the quarter stadium and extended B 
fields, Ref. ^6| gives k ~ 0.3. In this study we observed that k is a relatively robust number, 
since very different shapes of the conformal billiard give values of k that differ at most by 
60%. Therefore we believe that with help of Eq. (|3.9|) and taking k to be of order unity, one 
can estimate the magnitude of <f> c for other chaotic billiards. 



18 



IV. CONCLUSIONS AND DISCUSSION 



In this paper, we have proposed that the universal form of the parametric correlator of 
conductance peak heights indicates the chaotic nature of the electron dynamics in quantum 
dots in the Coulomb blockade regime. In experiments, the simplest parameter to vary is an 
external magnetic field. Whereas random matrix theory provides the universal form of the 
correlation function, the nonuniversal field scale can be understood in simple semiclassical 
terms: it is related to the average winding number per unit of time of periodic orbits 
bouncing between the confining walls of the quantum dot. This field scale is rather sensitive 
to the geometry of the dot and the Fermi energy. To compare our analytical and numerical 
predictions against the experimental result, the magnetic flux through the dot has to be 
larger than one quantum unit of flux h/e, but such that the cyclotron radius is much larger 
than the dot diameter. The former condition assures that time-reversal symmetry is broken; 
the latter implies that the bending of classical trajectories is mainly due to scattering by 
the boundaries. We point out that electron-electron correlations can be taken into account 
by assuming that the single-particle spectrum results from a self-consistent (Hartree-Fock) 
treatment.^ 

Our theoretical prediction for the correlator of peak heights is based on the hypothesis 
that the statistical properties of the system Hamiltonian can be described by random ma- 
trix theory. Although we could only derive analytically expressions for the limit of small 
field variations, the complete form of the correlation function was obtained by numerical 
simulations of large Gaussian matrices. We have compared the random matrix results with 
the exact correlator obtained from the conformal billiard after averaging over energy and 
shape deformation. The agreement found was good, given the limitations imposed by the 
size of the data set. In addition to that, we found that the result of the classical calculation 
for the magnitude of the field scale and its dependence on energy matches moderately well 
the quantum result. 

For experimental tests of our theory, it is important to note that dephasing in the small 
quantum dot has to be kept low enough, namely, the dephasing length has to be larger than 
the system size. Also, the parametric results presented for the GUE case were illustrated in a 
dynamical model in which only the magnetic field was varied. This was very convenient, since 
it allowed us to the relate the field scaling parameter to the underlying classical dynamics 
using a well-developed formalism. Alternatively, the external parameter could be taken as 
the shape variation. For B = 0, the conductance autocorrelation function should then follow 
our numerical results for the GOE case. Or, by varying the shape with B > B crossover , the 
conductance peak correlator should be given by our GUE results. (While finishing this 
work we learned that such an experiment has already been performed.^) Unfortunately, the 
semiclassical analysis for these situations is more difficult than the one presented here and 
still remains an open problem. 

Some billiard geometries have well-pronounced short periodic orbits which, for insufficient 
averaging over energy and magnetic field, can lead to strong nonuniversal features to the 
curves presented in this work. We believe that overcoming this problem will be one of 
the most serious challenges for the experiments. In particular, to verify experimentally the 
asymptotic behavior of c g (x) should be very difficult. Since c g (x) involves the subtraction 
of two numbers and for x ^> 1 these numbers become very close, statistical (nonuniversal) 
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fluctuations can easily drive the experimental c g (x) below zero. Another cause of deviations 
from the predicted universal behavior, presumably weaker, is the existence of wave function 
correlations which extend over the dot. In other words, our findings assume that the channels 
at different leads are decorrelated, which may not be completely true if the dot size is not 
much larger than the electron wavelength. The agreement of our random matrix results 
with the numerical calculations using the conformal billiard supports the assumption of 
independent channels for N > 200. Smaller systems should be more influenced by short or 
direct orbit effects. We should also mention that if there were strong correlations between 
heights of neighboring peaks of a given sequence, they would influence mostly the x ^> 1 
region of the correlator. This is because averaging over a finite range of magnetic field 
usually yields less statistics for large field differences and, consequently, more pronounced 
data correlation effects. 

Lastly, we point out the fact that, independently of previous considerations, another very 
interesting experiment is a direct measurement of c v {x). Despite the very large Coulomb 
energy, which makes the conductance peak spacing very regular at first sight very, a careful 
experiment should be able to observe the small fluctuations of the peak position as a function 
of an applied magnetic field. This will give direct information about the single-particle level 
dispersion or, equivalently, c v (x). It will also provide a direct test of our estimate of the flux 
correlation scale <p c . 

Note added. After the submission of this manuscript we learned of similar work by 
Alhassid et alM 
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APPENDIX A: EVALUATION OF (det(H 2 ) tx{H~ 2 )) H 



This Appendix is devoted to a rather detailed evaluation of the average over the deter- 
minant shown in Eq. (|2.20| ). Let us call 



c N = {det(H 2 )tr(H^)) H 



(Al) 



where H is now an iV x N GUE matrix. First, we notice that there is a more convenient 
way to express this average, namely, 



c N 



df{a) 
da 



(A2) 



a=0 
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where the generating function /(a) = (det(H 2 + aILjv))j? can be evaluated by the fermionic 
methodS 

f(a) = (det(H + at N ) det(H - al N )) H = (J d[x) exp[-x f (# ® U 2 + aljv ® > 

(A3) 

where a = —a 2 , L = diag(l, — 1), and x T — (Xi 1X2)1 with xi an d X2 representing N- 
component fermionic vectors. Averaging over the GUE matrix H we find that 



(exp(-x f # ® &2X))h = exp 
where u is the following 2x2 matrix 



A 2 



2N 



-trur 



(A4) 



u 



X1X1 X2X1 

X1X2 X2X2 



(A5) 



The quartic term can be decoupled by a Hubbard-Stratonovich transformation, namely, 



d[Q] exp 



-ytr(Q 2 ) - i\ti(Qu) 



2tt 



exp 



A 5 



2N 



-tr m z 



(A6) 



where Q is a 2 x 2 Hermitian matrix and <i[Q] = dQ\\dQ22dQvidQ2\- As a result, we have 



/(«) 



N\ 2 
2W 



exp 



-ftr(g 2 : 



d[x] exp[-x f (all7v ® £ + i\t N ® (A7) 



The Gaussian integral over the fermionic variables can be easily carried out, yielding 
f{a) = (£) 2 / ^P {-y tr(Q 2 ) + Mr[ln(aL + zAQ)]} . 



(A8) 



When iV ^> 1 the above integral over Q can be evaluated by the saddle-point approxima- 
tion (which becomes exact in the limit N — > 00). For this purpose we first separate angular 
and radial components of Q, namely, Q = T^qT, where q = diag(gi, g 2 ) and T is an SU(2) 
matrix. The differential breaks up into d[Q] = dfj,(T)J(q)d[q), where J(q) = 7r(gi — g 2 ) 2 is the 
Jacobian of the transformation, d[q] = dqidq 2 , and d/i(T) is the group measure normalized 
to unity. This yields 



N 
2^ 



d[q]J(q)exp(-NF[q]) / d/x(T) exp 



iaN 



tr(q- x TLT*) 



(A9) 



where F[q] = (l/2)tr(g 2 ) — tr[ln(zAg)] and we have only kept terms to lowest order in a. 
The saddle-point expansion now involves only the radial part of the action: 



F[q] = 1-2 ln(A) - ^{a x + a 2 ) + 5qf + 5q 2 2 + 0(5q 3 ), 



(A10) 
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where = cr 12 + 8qi,2 and a\ 2 = 1. The relevant saddle points correspond to <J\ = —02, 
resulting in 



f(a) = N (¥\ J dti(T) exp L l ^L t i(LTLT^ 



(All) 



The integral over the SU(2) manifold can be evaluated through the well-known Itzykson- 
Zuber formula,^ which in a simplified form reads 

J M T) eMMLTLT^ = ^ffl . (A12) 
with I12 denoting the eigenvalues of L. Hence, 

(the factor of 2 takes into account the double saddle point). Finally, we obtain 

C ~ = W 7 • (A14) 



APPENDIX B: THE ASYMPTOTIC LIMITS OF C r (X) 

The small-X asymptotics of the correlator of decay widths can be determined by the 
same method used in Sec. [11 A| for the conductance peak height correlator. Beginning with 
the definition presented in Eq. ([2.271) , we expand T ku (X ± X/2) up to first order in X [see 
Eq. ( |2.5|) ]. The zeroth-order term of Cr{X) is then given by 

Cr(0) = (r 2 ) - (U 2 = (r) 2 (Bl) 

for the unitary ensemble. An expression analogous to Eq. ( |2.8|) is used to write the second- 
order coefficient of Cr{X) in terms of the amplitudes A„, namely, 

Cf (0) = -^<A 2 >. (B2) 

Carrying out the average over the matrix elements of the external perturbation U [see 
Eq. flTT5D1, we find that 

C?(°) = -£(£^y (B3) 
We now average separately over the eigenvalues and partial widths and obtain 



27r 2 (n 2 cr 2 X 2 
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Upon rescaling both C r (X) -> c r (x) = C r (X)/(T} 2 and X -> x = X(a 2 /iVA 2 ), we arrive 
at 

Cr ( x ) = l-l^_ + 0(x 4 ). (B5) 
3 

The large- X asymptotics of Cr(X) can be inferred from the asymptotics of another 
correlator, namely, 

P(X,E) = n 2 \Mr;X 1 )\ 2 \Mr;X 2 )\ 2 6(E 1 - e fl (X l ))5(E 2 - s v (X 2 ))) - ij (B6) 



where Q is the system volume, X\p = X =F X/2, and £^2 = E E/2. Recall that the wave 
function intensities are proportional to the decay widths T^(Xi) and r„(X 2 ) for pointlike 
contacts. At E\ = E 2 and large X, the interlevel correlations are secondary to intralevel 
ones; as a result, the S function in Eq. ( |B6| ) acts as a Kronecker 5, causing P(X, 0) and 
Cr(X) to coincide (up to a prefactor equal to A 2 ) to leading order in 0{1/X). 

Let us for convenience assume a finite size space basis to represent the system Hamilto- 
nian. We can then reduce Eq. ( p6|) to 

N\ 2 „ . „ „ . „„ . 1 



P(X, E) = (- j (Im [G{E X + ze; X 1 )] kk Im [G(E 2 - ze; X 2 )] kk ) - — (B7) 

with G(E;X) = [E — H(X)]^ 1 and e — ► + . The above expression can be rewritten in the 
more convenient form 

P(X, E) = — ( Re ([GiE, + ie; X x )] kk [G(E 2 - ze; X 2 )} kk ) - -L . (B8) 

In general, an expression like Eq. (|B8| ) requires the evaluation of the following quantity: 

D klmn (E,X) = ([G(E 1 + ie;X 1 )] kl [G(E 2 -ie;X 2 )] mn ) . (B9) 

The correlator D k i mn {E, X) can be calculated exactly in the zero-mode approximation of 
the supersymmetric nonlinear a mo deli (or, equivalently, in the RMT framework). This 
calculation is standard nowadays (for a recent review, see Ref. |52|) and has already been 
presented in the literature.il Here we will only mention the resulting expression for the 
unitary ensemble, which is 

Dklmn(E,X) = (jy^J [Skl5 mn ~ 5 k l5 mn k(uJ, x) 

-5 kn 5i m n(u,x)) , (BIO) 

where 

dAij rfA 2 exp[2 7 ri(cj/2 + ^)(A 1 -A 2 )-(7rV/2)(A 2 -A 2 )J (Bll) 

and 



23 



n(u, x) = I dXt I 1 d\ 2 ( Xl+X * ] eX p \2m(uo/2 + ir))(Xi - A 2 ) - (ttV /2)(A? - A^) 



1 i-i \ Ai — A 



2 , 



(B12) 

When writing these equations we have rescaled the variables to E/A — u, e/A — rj, and 



NA^tY{U 2 )X/n 2 = x. Going back to Eq. ( pf) , we arrive at 

p{x, uj) = A 2 P(E, X) = ^Re [k{uj, x) + n(u, x)} . (B13) 

Since k(0, x) — > 2/ (irx) A and n(0, x) — > 2/(nx) 2 as x — ► oo, we have that p(x, 0) — > l/(7ra;) 2 
in the same limit. Therefore, we expect that 

cr(x) ™ -\ 2 . (B14) 

Finally, we remark that the x 3> 1 universal asymptotics of both k(0,x) and n(0,x) can 
also be obtained by the diagrammatic perturbation theory of disordered metals expressed 
in terms of diffuson modes. 



APPENDIX C: EVALUATION OF B(s) 



In this Appendix we give a schematic description of the calculation of B(s). The starting 
point is Eq. ( |2.37| ). Here we go through the same steps of Sec. |11 A| to evaluate A £ [see 
Eq. ( p.l8| )1. First we fix the reference eigenvalue e v to the center of the spectrum, obtaining 



B(k) 



I- 



N 



(kA/2wf 



-1\ 



(CI) 



{e v } 



Next, we rephrase this expression in terms of an average over a spectral determinant, namely, 

; det (# 4 ) 



B(k) = oat 



,det [H 2 + (k A/2 




(C2) 



where is a constant [such that B(0) = 1] and the average is performed over a 
(N — 1) x (iV — 1) GUE matrix H. The appearance of determinants in both numerator 
and denominator in Eq. (|C2j ) makes its evaluation technically more difficult than Eq. ( |A1| ). 
It is necessary to introduce not only four anticommuting auxiliary variables, but also two 
commuting (complex) ones. The resulting symmetry group is U(l, 1|4) (the pseudo-unitarity 
is due to the structure of the denominator). Fortunately, a general solution for such graded 
symmetry problems has been recently worked out.c3 The derivation is a nontrivial general- 
ization of the method of Appendix 0. For an expression with the structure of Eq. (|C2| ), one 
arrives at the following formula:^ 
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' II- = idet(ff-m 3 -A) \ 
^det(H - iaA) det(H + iaA) / H 



Ae 



—2wa 



a 



£ 

{™j} 



(ia — 777.3) {iot — mi) 
(m 3 - mi)(m 3 - m 2 ) 



x (m + mi)(m + m 2 ) c ^ (mi+m2 - m3 - m4 ) 
(7774 — 777i)(m4 — m 2 ) 



(C3) 



where the sum runs over all six nonequivalent combinations of pairs of rrij, a > 0, and A is 
an unspecified constant. To get B(k) we need to take the limit rrij — > for all j = 1, 2, 3, 4 
at a given order. After some algebra, one finds that 



B(k) = e — (24 + 24A; + 8A; 2 + A; 3 ) . 
Finally, inverse Fourier transforming the above expression, we arrive at 



(C4) 



dk 

35 + 14s 2 
12tt(1 + 



tfcs 



3s 4 



(C5) 
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FIGURES 
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FIG. 1. The rescaled correlator of conductance peak heights (O) obtained from the Hermitian 
random matrix simulations (unitary ensemble). For comparison, the inset also shows the large x 
asymptotics in the correlator of decay widths (A). The solid lines are the analytical predictions 
for the asymptotic behaviors (when known) and the dashed line corresponds to the fitted curve 
c g (x) = 0.735(7rx) -2 . Statistical error bars are too small to be seen. The arrow indicates the 
correlation width at half maximum height. 
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FIG. 2. The rescaled correlators of conductance peak heights and decay widths from real ran- 
dom matrix simulations (orthogonal ensemble), following the same conventions as Fig. [l]. The 
dashed line corresponds to the fitting c g (x) = 1.214(7rx)~ 2 and the solid line is the theoretical 
prediction for the i»l asymptotics of cr(x). 
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FIG. 3. The Fourier transform of the correlator of peak heights for the GUE (O) and GOE (A) 
ensembles. The dashed line is the curve f(k) = 0.125e~ fe / 2 ' 75 representing the Fourier transform 
of a Lorentzian fitted to the GUE data. 
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FIG. 5. The decay width correlation function for the conformal billiard (O) for a fixed shape 
= c = 0.2 and S = 7r/3) and averaged over 76 values of flux and 200 energy levels. The full line 
the result of the GUE simulation. 
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FIG. 6. The peak height correlation function for the conformal billiard (O) f° r a fixed shape 
(6 = c = 0.2 and S = 7r/3) averaged over 76 values of flux and 200 energy levels. The full line is 
the result of the GUE simulation. 
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FIG. 7. Numerical estimate of the winding number variance w 2 as a function of orbit time T 
in the classical billiard with unit area, vf = 1, and 5 = tt/3. 
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FIG. 8. The distribution of winding numbers (solid line) for T = 150 in the classical billiard 
with unit area, vp = 1, and 5 = tt/3. The dotted line is a Gaussian curve with variance given by 
Eq. | 
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FIG. 9. The numerical coefficient k as a function of geometry for the classical conformal billiard. 
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FIG. 10. C„(0) as a function ofyN for the conformal billiard, with iV as the eigenstate number. 
The symbols indicate the data obtained from the exact numerical diagonalization (fully quantum) 
for different geometries. The solid line is the total average over all data: C v (0) = 1.202^/N. The 
dashed and dotted lines are the semiclassical estimates [Eq. (^)] for 5 = ir/2 (k = 0.253) and 
5 = 7r/3 (k = 0.314), respectively. 



37 



